library(tidyverse)
library(texreg)
library(lmtest)
library(sandwich)
library(ggplot2)
library(haven)

#### Table 1: Regression Estimates of Support for War (PureSpectrum Survey) ####

# load the cleaned dataset 
df_clean <- readRDS("PureSpectrum/clean_data_PureSpectrum.rds")

# run the regression models
m1 <- lm(attack ~ hmrts + alliance, data = df_clean)
m2 <- lm(attack ~ alliance * hmrts, data = df_clean)
m3 <- lm(attack ~ alliance*hmrts + 
                male + age_cat + edu4 + inc_10k, data = df_clean)

texreg(l = list(m1, m2, m3),
       reorder.coef= c(2, 3, 4, 5, 6, 7, 8, 1),
       custom.coef.names = c("(Intercept)", "Violating Human Rights", "U.S. Military Alliance",
                             "Violating Human Rights $\\times$ U.S. Military Alliance",
                             "Male", "Age", "Education", "Income"),
       stars = c(0.01, 0.05, 0.1),
       digits = 2,
       caption = "Regression Estimates of Support for War (PureSpectrum Survey)",
       caption.above = T,
       include.ci = F,
       include.rmse =F,
       include.rsq = F,
       include.adjrs = F,
       custom.note = "",
       label = "tab:main",
       fontsize = "small") %>%
  gsub(".begin.center.", "\\\\centering", .) %>%
  gsub(".end.center.", "", .)


#### Figure 1: Impact of Treatments on Support for War (95% Confidence Intervals) ####

### Plot (1a) Support for Attack
# calculate the mean support for attack in each treatment group
df_ate <- df_clean %>% 
  group_by(exp_4, alliance, hmrts) %>%
  summarise(ate = mean(attack, na.rm = TRUE),
            n = n(),
            sd = sd(attack, na.rm = TRUE),
            se = sd(attack, na.rm = TRUE) / sqrt(n)) %>%
  mutate(ci_low = ate - 1.96*se,
         ci_high = ate + 1.96*se)

# plot the mean support for attack in each treatment group and the 95% CI
ggplot(df_ate, aes(x = factor(exp_4, level=c('1', '3', '2', '4')), y = ate,
                   color=factor(exp_4))) + 
  theme_classic() +
  geom_point()+
  geom_errorbar(aes(ymin=ci_low, ymax=ci_high), width=.1,
                position=position_dodge(0.05)) +
  scale_x_discrete(labels= c('Non-Ally,Violates HR', 'Ally,Violates HR',
                             'Non-Ally,Protects HR', "Ally,Protects HR")) +
  scale_color_manual(values=c('#999999', 'black', 'black', 'black')) +
  theme(legend.position = "none") +
  labs(x = "", y = "Mean Support for Attack \n (% point)", size = 12) +
  geom_text(aes(label=round(ate, 1)), position=position_dodge(width=0.9), 
            vjust=.5, hjust = -.35) +
  theme(axis.text.x = element_text(angle = 20, hjust = 0.5, vjust = 0.5, size = 12),
        axis.title.y = element_text(size=12))

# save the plot
ggsave("ate-1.pdf", width = 6, height = 4)


### Plot (1b) Average Treatment Effect on Support for Attack
# calculate the difference in support for attack between 2-4 against baseline condition
est <- rep(NA, 4)
ci_low <- rep(NA, 4)
ci_high <- rep(NA, 4)
se <- rep(NA, 4)

for(i in 2:4){
  test <- t.test(df_clean$attack[df_clean$exp_4==i], 
                 df_clean$attack[df_clean$exp_4==1])
  est[i] <- test[["estimate"]][["mean of x"]] - test[["estimate"]][["mean of y"]]
  ci_low[i] <- test[["conf.int"]][1]
  ci_high[i] <- test[["conf.int"]][2]
  se[i] <- test[["stderr"]]
}

df_ate_diff <- data.frame(exp_4 = df_ate$exp_4, est, ci_low, ci_high, se)
df_ate_diff <- df_ate_diff[-1, ]

# plot the differences and 95% CI
ggplot(df_ate_diff, aes(x = factor(exp_4, level=c('3', '2', '4')), y = est)) + 
  theme_classic() +
  geom_point()+
  geom_errorbar(aes(ymin=ci_low, ymax=ci_high), width=.1,
                position=position_dodge(0.05)) +
  scale_x_discrete(labels= c('Ally,Violates HR', 'Non-Ally,Protects HR', 
                             'Ally,Protects HR')) +
  labs(x = "", y = "Average Treatment Effect \n (% point)") +
  geom_text(aes(label=round(est, 1)), position=position_dodge(width=0.9), 
            vjust=.5, hjust = -.35) +
  geom_hline(yintercept = 0, linetype="dotted") +
  theme(axis.text.x = element_text(angle = 20, hjust = 0.5, vjust = 0.5, size = 12),
        axis.title.y = element_text(size=12))

ggsave("ate-diff-1.pdf", width = 6, height = 4)

